# Dickstein, Ho, and Mark (2023)
# Estimating demand models with randomly drawn households

# Preliminaries
setwd("../../../../library")
source("PreliminariesCode.R")
library(parallel)
set.seed(1)

#Establish seed setting for parralel processing: 
RNGkind("L'Ecuyer-CMRG")
generate_randomseed <- function(run){
  s <<- .Random.seed
  if(run > 1){
    for (i in 1:(run - 1)) {
      s <<- nextRNGStream(s)
    }
  }
  .Random.seed <<- s
}

## Loading cleaning function: 
source(paste0(project_folder, "/library/DA01QuanModifyExplodedDFunc.R"))

# Setting options
ravec <- 1:7
yrvec <- 2014:2016
specs_location <- paste0(project_folder, "/analysis/estimationcode/specs/lowrisk_simplemh_censor0.2_switcher")
specs_name <- "lowrisk_simplemh_censor0.2_switcher"
testing <- T
fixed_cost_censor <- 0.2
sumrisk_cutoff <- 70
est_type <- "switcher_only"
check_derivatives <- FALSE

# Loading data
explodeddata <- fread(paste0(project_folder, "/data/explodeddata.csv"))
cgquinsvec <- as.numeric(scan(file = paste0(project_folder, "/data/orig/acg_positive_quintiles"), sep = ",", what = "character")[2:7])

# Making the cleaning easier by keeping only switcher and small other sample: 
explodeddata_switcher <- explodeddata[switcher == 1]
explodeddata_balancer <- explodeddata[age == 61]

# Adding the function that takes the explodeddata set, re-samples, then re-estimates demand
bootstrap_estimate <- function(k){
  print(paste0("Starting bootstrap draw", k))
  
  ## *BOOTSTRAP* creating a random sample with replacement
  samp <- sample(unique(explodeddata_switcher[switcher == 1]$hhid), length(unique(explodeddata_switcher[switcher == 1]$hhid)), replace = TRUE)
  setkey(explodeddata_switcher, "hhid")
  explodeddata_k <- rbind(explodeddata_switcher[J(samp), allow.cartesian = TRUE], explodeddata_balancer)

  ## Cleaning
  print(paste0("Modifying... ravec: ", paste(ravec, collapse = " and "), "yrvec: ", paste(yrvec, collapse = " and ")))
  data <- mod_exploded_data(
    explodeddata_k,   
    paste0(project_folder, "/data"),
    highmeanthresh = 0.59616,
    lowsumthresh = 0.2575,
    highsumthresh = 1.76329,
    highmaxthresh = 3.37157,
    medincome = 2.5,
    acgquins = acgquinsvec,
    winsorbound = 89.3707340916662,
    ravec = ravec,
    yrvec = yrvec)
  
  ## Change where to save the outcome
  output_path <- paste0(specs_location, "/bootstrap/output_", k)
  dir.create(output_path, recursive = T)
  print(paste0("output path:", output_path))
  
  # Restrict data to switcher sample
  data <- data[switcher == 1]
  
  # Set/determine variables
  not_zero <- function(x) ifelse(is.numeric(x), sum(x, na.rm = T) != 0, F)
  nums <- unlist(lapply(data, not_zero))
  payervars <- grep("PAYER_ID_", names(data[, ..nums]), value=TRUE)
  var_names <- list(
    W1_names = c("low_sum_risk", "sum_concurrent_risk",
                 "withkids", "over50"),
    W2_names = c("cons"),
    W3_names = c("cons"),
    Z_names  = c(payervars, "silver_costshare")
  )
  var_switch <- list(
    W1_names = rep(0, length(var_names$W1_names)),
    W2_names = rep(0, length(var_names$W2_names)),
    W3_names = rep(0, length(var_names$W3_names)),
    Z_names = rep(0, length(var_names$Z_names))
  )
  
  # Calling estimate demand!
  source(paste0(specs_location, "/../../estimatedemand.R"), local = TRUE)
  rm(data)
}  

# START LOOP

# How many cores do you have available?
numCores <- detectCores()
k_vec <- 1:20
# Non-parallel processing version:
for(k0 in k_vec){
  set.seed(1)
  generate_randomseed(k0)
  bootstrap_estimate(k0)
}
# Parallel processing version:
#mclapply(k_vec, bootstrap_estimate, mc.cores = numCores - 1, mc.set.seed = T)

# END LOOP


